{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5678b63a-3352-4256-8691-d9578a292fa9",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Estimate structural parameters: sensitivity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fe55b8d-1499-4260-bbe1-1eba710baf7f",
   "metadata": {
    "tags": []
   },
   "source": [
    "## The baseline setting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "56a0e031-ed12-4bb1-9820-7183d7a4c0e4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "from fn_simulate import *\n",
    "\n",
    "np.set_printoptions(precision=4)\n",
    "\n",
    "# randomization variables for simulation\n",
    "#import secrets\n",
    "#secrets.randbits(128) \n",
    "seed = 84708628852577542122415832887111152745\n",
    "\n",
    "rng = np.random.default_rng(seed)\n",
    "\n",
    "n_draws = 1000\n",
    "\n",
    "param_random = {'n_draws' : n_draws, \n",
    "                'seed' : seed}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ca330fd5-a3cf-4d09-a8af-b0820508d2f9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# calibrated parameters: baseline values\n",
    "\n",
    "    # income\n",
    "mu_lny = 9.713\n",
    "sig_lny = 0.497\n",
    "\n",
    "mu_eps = 0\n",
    "\n",
    "    # prices\n",
    "pi_s = 1\n",
    "pi_q = 0\n",
    "tau = 0.075\n",
    "\n",
    "gamma = 0.424\n",
    "\n",
    "param_given = {'mu_lny' : mu_lny, \n",
    "               'sig_lny' : sig_lny,\n",
    "               'mu_eps': mu_eps,\n",
    "               'pi_s' : pi_s, \n",
    "               'pi_q' : pi_q, \n",
    "               'tau' : tau,\n",
    "               'gamma': gamma}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bf3cbf04-eda7-4325-a452-167a20269d6c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# estimated parameters\n",
    "\n",
    "# initial:\n",
    "pi_n = 1.467219\n",
    "pi_nq = 0.110516\n",
    "\n",
    "theta = 0.927190\n",
    "alpha = 0.611335\n",
    "rho = -7.924372\n",
    "\n",
    "sig_eps = 0.299971\n",
    "eps_cutoff = -0.139612\n",
    "\n",
    "cor_lny_eps = -0.010366\n",
    "\n",
    "\n",
    "param_choice = np.array( \n",
    "    [pi_n, pi_nq, \n",
    "     theta, alpha, rho,\n",
    "     sig_eps, eps_cutoff, cor_lny_eps] \n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bccf8ec9-4482-4e72-8528-7f09b1a5f5e3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# moments in the real data\n",
    "\n",
    "# moments observed in real data\n",
    "beta_A_2to3_real = -0.157  \n",
    "beta_B_2to3_real = 0.412 \n",
    "\n",
    "PA_2to3_real = 0.212\n",
    "PB_2to3_real = 0.121\n",
    "\n",
    "exp_n_AB_real = 0.214\n",
    "exp_nq_AB_real = 0.111\n",
    "\n",
    "cor_lny_n3_real = 0.128\n",
    "cor_lny_lnq_AB_real = 0.296\n",
    "\n",
    "moments_real = np.array( \n",
    "    [beta_A_2to3_real, beta_B_2to3_real, PA_2to3_real, PB_2to3_real, \n",
    "     exp_n_AB_real, exp_nq_AB_real, cor_lny_n3_real, cor_lny_lnq_AB_real] \n",
    ")\n",
    "\n",
    "# weights: 1/moments_real^2 as in Baudin, de la Croix, & Gobbi (2015, AER)\n",
    "W = np.diag( 1 / np.square(moments_real) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8c43c491-3281-4709-ac76-a44f940e39f7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Display for each simulation: param_choice and value of the objective function\n",
    "disp = True\n",
    "\n",
    "# Define the objective function to be minimized\n",
    "def objective_function(param_choice): # choice variables should be a numpy array\n",
    "    \n",
    "    return simulate_objective(simulate_moments, param_choice, param_given, param_random, moments_real, W, disp = disp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d3e7e367-fd4c-421f-baf3-83e7519d20f1",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------------------------------------------------------\n",
      "Evaluating the objective function via a simulation ... \n",
      "Parameter names: pi_n pi_nq theta alpha rho sig_eps eps_cutoff cor_lny_eps\n",
      "param_choice: [ 1.4672  0.1105  0.9272  0.6113 -7.9244  0.3    -0.1396 -0.0104]\n",
      "Real moments: [-0.157  0.412  0.212  0.121  0.214  0.111  0.128  0.296]\n",
      "Simu moments: [-0.157   0.412   0.2121  0.1209  0.2127  0.111   0.1286  0.296 ]\n",
      "Distance    : [ 3.1756e-05 -2.5513e-05  8.9077e-05 -1.0923e-04 -1.2806e-03  3.2279e-05\n",
      "  6.4374e-04  4.0237e-05]\n",
      "Objective function value:  6.224028519632813e-05\n",
      "Simulation time: 4.09524393081665 seconds\n",
      "-------------------------------------------------------\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAAUCAYAAABmtzNiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAABJ0AAASdAHeZh94AAAJU0lEQVR4nO2cf7BVVRXHPw9RIVT8FTJpkZCEiomO4o9MQYgSzXykM44jJU4ak0VPZUwd9cuqMXEUf9tU2ACaU6Ap+QPxR5CkKMyohUgJAQ/LEvPnoIKK2h9rH99+h3PuPffc+8CG+525s985e6+9115nnb3XWnud1/LRRx/RRBNNbN3ovqUZaKKJJroWZjYJUOr2Wkl9k4vmQtBEE1sHngeGRdcfxJXNhaCJJrYObJT0Ul5l5kJgZl8B2oAjgV2B14BngeskzcnrzMx2A1qB44EDgD2B9wLtNGCapA/rpckZeyxwa7g8S9ItOe32An4CfB3YDfgPMBswSa83gi8zOx74EbBfNMZTwDWSnsho3w70y5laJxMuojkZOAYYAhwI7AjcLun0HJ5agHHA2cD+wDb4LjENuFnSB1l0gbawPpjZlcAhwEBgd2A9sAaX8U2SXo3abk4Z1ySvelBmrKJ6WQf6m9mLuHyfBC6S1J5Udstg6BJgAXA0MBeYAtwL7EJn0yILpwBTgcOARcB1wO+BwcAtwKygkPXSpHn+LHAj8FaVdgNwZRkHLAauBVbhCvVEUMy6+Aovwn3Awbj8rgeeBr4JPG5mecrwJmAZv6tz2l8C/ABXthcrzTtgBvBrYG9gZpjbdoG/mXkyLqEP5wK9gIdD37cDG4FJwJLwrBJsThnXKq96UNNYNeplGSwCxuKLzFlAn9DvrkmDThaBmZ0C/BR4BBgjaV2qftsqAy4HTgTuj1dyM7sYn+C3gDH4w66HJuapBd89XgXuAiZW4O/nuBAmSLox6uMaXIEvB8aX5cvM+obx1wJfkvRyVDccmIev+r/J4O0NSZMq8J7GucC/gH/gu8/8vIZmdhKuCKuBoZJeCfe3BWaFuXwHmJ6iK6MPO0nakMHD5cDFwEXA98PtzSnjwvJqAGodqxa9zAv+pTFc0p8AJD0Q3X/WzJ7AF5ozgGsgsgjMrBtwJfAOcFr6oYcO3680sqR5ku5Nm3PBN/lFuBxWL00KE4Bj8dX07bxGZtYfGAW0AzenWQ+0Y82sVx189cNluihW0EA3H1gHfLrCXApD0nxJKyQVOf8dE8opySIQ+ngfuDRc/jAmKKsPWYtAwKxQ7hO13WwyrlFedaGWsWrVy4CbgH2r/BZX4O9tYBnRs4gtgiNxs/FO4PXggw0GNgCLs/yuGpEozcZG0ZjZvsBk4HpJC8zs2Ap9JXUPZSjeOjN7HH8ghwN/LMnXCtwHG2pmu8cvnZkdjfuKs3P63D6YtJ/DH/4SYEEl370GJDGGVRl1yb2DzWxnSW+E60brwzdCuaRg+66Q8ScRNetlmPMrlISZ9QAGEVkq8UJwaCjX4v7WASniBcDJkv5bYuDuwLfD5dxG0IT624AXcJOzGr4YyuU59StwgQ+kwkJQiS9Jr5nZj3Fza5mZzcZdlgG4Cfww8L2crvvi84mx2szGSXo0j5+CSJRm74y6/tHfg/BAEtSpD2Y2EdgB6I0HD4/CF4HJ1ZjtQhmXhpnth8dUnpe0voFdN0QvK8HMrsbjOi/gLsileBxnRtImDhb2CeV4oCcwEl9dBwMP4sGiO8owgj/8wcAcSQ82iOYy4CDgjIIPpnco38ypT+7vXA9fkq7DTfHueGDmQjwo9k9getqcDZgGjMAXg174S/dL4PPAA2Z2YBWequG+UJ4XB4jCC2dRu12iv+vVh4m4aduGLwJzgVEFN5KukHG9mAM8Q8eL2yg0Si8rYS/gt/gp0V3Au8DhktYkDWKLYJtQtuAr/V/D9XNm1oqvWMeY2RG1mIVmNgE4H/g7HrCqm8bMhuJWwJQGuCwJkuh0rl9XZC5mdgHwM+AG3Jd7Cd9prwBuN7Mhki6IaSRZqpulwHgzeyuMNwk/ZiuL3wGnA8fhu+g9uO8/Et9JV+D+YuyG1KUPyZGnme2BuxmTgWfM7ARJT+cx2lUy/j9GVb2sBkmnVmsTWwTJWeWq6KEnHa3HdwGAoUUZMLNz8KOdZXgU87V6aSKXYDkdga4iSFbW3jn1O6Xa1cRXaDMMD7DdI+k8SaskvRMUvxU/Sjo/BIiKIAmYHV2wfSaC73kivku/hL9gZ+KR7aNw0xog3kkbog+S1kq6Gzdvd6Mj12MTbCEZF8WZeJwjK85SD+rSy0YhtgieD+UbOW0TxehZpGMza8PPQ5cCI4qYawVpdsD9JYANZunNFICpZjYVDyK2hXvJ/AZmEdARQd3EV6thLieEcpPjIknvmNliXFkPophCJeP0qtiqACRtxHMApsT3zawnft69HnguqmqoPkhaY2bLgCHpIF/go40tI+NCkDSvUX2lUFovG4nYIliAR2j3MbPtMtoODmV7tU5DMOda4C/4yl5kEShK8y6eGJP1eya0eSxcxyZrojijwtFYPPaOwJfxl+HJVF0tc9k+lHlHhMn99yr0EeOIUDZ6F4oxFugBzEodBzZMHyJ8JpSdTkK2sIy3NErpZaPx8cBhhZ6JmyiXpRj6KvA13DyZG90fYGaD4sQSM7sU9wefwlf2qscctdBIWi/pu1k/4J7QbEa4NzOiWwk8hAfgzkmzgO+6t4Yz1rJz+XMozzazPVNzPA5/qBuAhdH9/eMAXnS/H+7/QnYCUk0ws50y7h2Kz+8tPAnnY5TUh0Eh4Sc9TreQUNQHWKjOqdxdLuNPMsroZVegJf5/BGbWB3gc+AIu8MV4AkcrHqw4TdIdUfv2UL+3pHYzS7LTPsBTfrP8mnZJ06M+aqbJQ5RxlfmtgXkq50JcIf8A/A1Pbx2Om15HKuTCl5xLN9x3HoknttyN++T74iZtC9Am6foUzxfiO8PqQDcAz7/vgUerWyV12uFCtuBJ4bIv/mKuouNFeUXSxKj9InxnWRrG2B8YjVtYY7Ki8yX0oQ24CrcmVuKxhz3w7Lr+QRYjJC0L7TeLjMvIqx6UeDaF9bKr0CnFWNLLZnYYnivdiicxrAPuB66QVM08Sc6pt8GPjbLwKJ1TWcvQlIKklWZ2CB0fd4zGP+64Af+4Iw5O1cyXpA/NbDS+sp+Ky/BT+Ec6c4AbJD2U6mM+fiR1EO4K9ML98sfwoOhtORlqQ/C04Bj96cgLWEPndOs7A0+n4379v/F8/smKPj6JUUIfHgF+he/KB+JHXm/jynxbmP+WkDHULq96UNNYNepll6Cl+R+KmmiiiU2+PmyiiSa2PvwPktxYOs4/OBoAAAAASUVORK5CYII=",
      "text/latex": [
       "$\\displaystyle 6.22402851963281 \\cdot 10^{-5}$"
      ],
      "text/plain": [
       "6.224028519632813e-05"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# a trial simulation to check execution time \n",
    "objective_function(param_choice)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f03a2f18-c570-41e9-8cc4-ae0204600dca",
   "metadata": {},
   "source": [
    "## Sensitivity to $\\tau$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "28071296-a0be-41cc-b598-88df8a03af73",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# numpy display precision\n",
    "np.set_printoptions(precision=4)\n",
    "\n",
    "# Display option param_choice and value of the objective function, for every simulation inside an estimation \n",
    "disp = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a4f46e94-d5a5-4d79-8087-1c8dc6c0a209",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****** The value of tau ******:  0.05\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.032431\n",
      "         Iterations: 415\n",
      "         Function evaluations: 731\n",
      "--- 3227.5089762210846 seconds ---\n",
      "****** The value of tau ******:  0.06\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.008640\n",
      "         Iterations: 380\n",
      "         Function evaluations: 715\n",
      "--- 4550.495681762695 seconds ---\n",
      "****** The value of tau ******:  0.07\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.000863\n",
      "         Iterations: 333\n",
      "         Function evaluations: 558\n",
      "--- 5382.171635389328 seconds ---\n",
      "****** The value of tau ******:  0.08\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.000769\n",
      "         Iterations: 260\n",
      "         Function evaluations: 484\n",
      "--- 3242.54869055748 seconds ---\n",
      "****** The value of tau ******:  0.09\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.005697\n",
      "         Iterations: 322\n",
      "         Function evaluations: 562\n",
      "--- 3628.715814113617 seconds ---\n",
      "****** The value of tau ******:  0.1\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.019256\n",
      "         Iterations: 311\n",
      "         Function evaluations: 553\n",
      "--- 3438.9420070648193 seconds ---\n"
     ]
    }
   ],
   "source": [
    "# list of tau for sensitivity tests\n",
    "tau_sens_list = [0.05, 0.06, 0.07, 0.08, 0.09, 0.10]\n",
    "\n",
    "# a container for optimized parameters\n",
    "tau_sens_param_choice_list = []\n",
    "\n",
    "# iterate over different tau\n",
    "for tau in tau_sens_list:\n",
    "    \n",
    "    print(\"****** The value of tau ******: \", tau)\n",
    "    start_time = time.time()\n",
    "    \n",
    "    # reset the value of tau in param_given, which will be passed to the optimization function\n",
    "    param_given['tau'] = tau\n",
    "    \n",
    "    # Define the objective function to be minimized\n",
    "    def objective_function(param_choice): # choice variables should be a numpy array\n",
    "\n",
    "        return simulate_objective(simulate_moments, param_choice, param_given, param_random, moments_real, W, disp = disp)\n",
    "\n",
    "    # Perform the optimization\n",
    "    result = minimize(objective_function, param_choice, method = 'Nelder-Mead', \n",
    "                      options = {'disp':True, 'fatol':0.0001, 'maxiter':1000}\n",
    "                     )\n",
    "\n",
    "    # store results\n",
    "    tau_sens_param_choice_list.append(result.x)\n",
    "    \n",
    "    print(\"--- %s seconds ---\" % (time.time() - start_time))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57432149-3b90-4ad2-b59b-fe892444d23b",
   "metadata": {},
   "source": [
    "## Sensitivity to $\\gamma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ffe9964a-d31b-42de-9c9e-c87a23405e5a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****** The value of gamma ******:  0.3\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.248166\n",
      "         Iterations: 672\n",
      "         Function evaluations: 1127\n",
      "--- 8129.228497982025 seconds ---\n",
      "****** The value of gamma ******:  0.35\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.082428\n",
      "         Iterations: 399\n",
      "         Function evaluations: 706\n",
      "--- 4788.962294816971 seconds ---\n",
      "****** The value of gamma ******:  0.4\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.032551\n",
      "         Iterations: 296\n",
      "         Function evaluations: 541\n",
      "--- 3575.5010941028595 seconds ---\n",
      "****** The value of gamma ******:  0.45\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.014556\n",
      "         Iterations: 299\n",
      "         Function evaluations: 539\n",
      "--- 3161.4071719646454 seconds ---\n",
      "****** The value of gamma ******:  0.5\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049970\n",
      "         Iterations: 209\n",
      "         Function evaluations: 409\n",
      "--- 2229.0562369823456 seconds ---\n",
      "****** The value of gamma ******:  0.55\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.066220\n",
      "         Iterations: 242\n",
      "         Function evaluations: 440\n",
      "--- 2330.8141322135925 seconds ---\n"
     ]
    }
   ],
   "source": [
    "# list of tau for sensitivity tests\n",
    "gamma_sens_list = [0.30, 0.35, 0.40, 0.45, 0.50, 0.55]\n",
    "\n",
    "# a container for optimized parameters\n",
    "gamma_sens_param_choice_list = []\n",
    "\n",
    "# iterate over different tau\n",
    "for gamma in gamma_sens_list:\n",
    "    \n",
    "    print(\"****** The value of gamma ******: \", gamma)\n",
    "    start_time = time.time()\n",
    "    \n",
    "    # reset the value of tau in param_given, which will be passed to the optimization function\n",
    "    param_given['gamma'] = gamma\n",
    "    \n",
    "    # Define the objective function to be minimized\n",
    "    def objective_function(param_choice): # choice variables should be a numpy array\n",
    "\n",
    "        return simulate_objective(simulate_moments, param_choice, param_given, param_random, moments_real, W, disp = disp)\n",
    "\n",
    "    # Perform the optimization\n",
    "    result = minimize(objective_function, param_choice, method = 'Nelder-Mead', \n",
    "                      options = {'disp':True, 'fatol':0.0001, 'maxiter':1000}\n",
    "                     )\n",
    "\n",
    "    # store results\n",
    "    gamma_sens_param_choice_list.append(result.x)\n",
    "    \n",
    "    print(\"--- %s seconds ---\" % (time.time() - start_time))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e20d9c4-c476-47c7-b869-f26aa9a36df6",
   "metadata": {},
   "source": [
    "## Store estimated parameters in the sensitivity tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d11abdef-84e9-4141-accf-5aa089839e3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_names = ['pi_n', 'pi_nq', 'theta', 'alpha', 'rho', 'sig_eps', 'eps_cutoff', 'cor_lny_eps']\n",
    "\n",
    "df_tau_sens_param_choice = pd.DataFrame(tau_sens_param_choice_list, columns = param_names)\n",
    "\n",
    "df_gamma_sens_param_choice = pd.DataFrame(gamma_sens_param_choice_list, columns = param_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b9da3e34-7236-4a45-9f0c-aa0f3cd03427",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_tau_sens_param_choice.to_csv('./df_tau_sens_param_choice.csv')\n",
    "\n",
    "df_gamma_sens_param_choice.to_csv('./df_gamma_sens_param_choice.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "929e3bd3-5768-4f30-b963-8422855198e8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:base] *",
   "language": "python",
   "name": "conda-base-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
